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Condiciones de estabilidad de la Optimización 
Dinámica Estocástica aplicada al cálculo del valor 
del agua de un embalse 


R. Chaer y P. Monzón 


Abstract—El presente trabajo analiza condiciones necesarias y 
suficientes para la estabilidad de un algoritmo de optimización 
dinámica estocástica para el cálculo del costo del agua en un 
sistema eléctrico con una planta de generación hidroeléctrica. 
Se estima la región de estabilidad y se presentan diferentes 
simulaciones para parámetros dentro y fuera de la misma. Se 
establece una relación entre el paso de integración temporal, la 
discretización espacial y los máximos flujos de entrada y salida 
del embalse. A efectos de mostrar ventajas y desventajas del 
método propuesto, se presentan simulaciones del sistema eléctrico 
uruguayo con cuatro centrales de generación hidroeléctrica. En 
este caso, las condiciones de estabilidad imponen un paso de 
integración temporal muy pequeño, lo que resulta en un tiempo 
de simulación muy grande. 

Index terms: Optimización estocástica, Despacho de potencia. 


I. INTRODUCCIÓN 


A optimización de la operación de un sistema hidrotérmico 

resulta ser un problema complejo. Esto se debe a la presencia 
de embalses y su consecuente operación, que involucra decisiones 
respecto a cuándo usar el agua y cuánto usar. Estas decisiones afectan 
no solo la operación actual, sino también la futura. 

El problema puede ser formulado como la minimización de una 
función objetivo: la suma del costo del combustible en las centrales 
térmicas y del costo de falla por no poder satisfacer la demanda 
de potencia y energía, con la restricción que imponen la red de 
transmisión y los límites operativos. Para su resolución, se divide 
el tiempo en sucesivos pasos o etapas. En cada etapa, se suponen 
conocidos los costos de producción de cada unidad térmica y y el 
costo de falla del sistema. Como el costo del agua en los embalses 
no tiene un costo explícito asignado, no es posible manejar un costo 
de producción asociado a las centrales hidroeléctricas. Sin embargo, 
hay que tener en cuenta que si el agua en el embalse se usa hoy, 
la disponibilidad futura puede verse comprometida. Por otro lado, 
preservar el agua hoy, puede redundar en menos costos en el futuro, 
ya que puede no ser necesario prender más centrales térmicas. El 
problema es tener elementos de base para una política de uso de los 
recursos acumulados que resulte en un compromiso favorable entre 
los costos actuales y futuros. 

En ese sentido, se enfrenta un problema de optimización: mini- 
mizar una función de costo sujeta a un conjunto de restricciones. Si 
bien hay varios métodos para atacar el problema, los más clásicos son 
el denominado Stochastic Dynamic Programming (SDP) y el Stochas- 
tic Dual Dynamic Programming (SDDP), que incorpora el problema 
dual. El SDP calcula la función de costo desde el futuro hacia atrás, 
hasta el momento actual. Para realizar este calculo, se realiza una 
doble discretización, tanto en el tiempo como en el espacio, para 
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todas las variables de estado del sistema. Esto lleva naturalmente a 
trabajar con gran número de variables y a la denominada maldición 
de la dimensionalidad de Bellman, por lo que este método no es 
aplicable para sistemas grandes [1], [2]. Por otro lado, el SDDP 
ataca el problema de la dimensionalidad a través de la resolución 
del problema dual realizando los denominados cortes de Benders 
para aproximar la función de costo. Una buena explicación de este 
procedimiento puede encontrarse en [3]. El método funciona bien 
con funciones convexas, pero esto no permite incorporar restricciones 
naturales, como los mínimos operativos de las centrales térmicas. La 
no convexidad se traduce en un gap que suele ser despreciable en 
sistemas grandes [4], pero puede ser importante en sistemas medianos 
y pequeños, en los que, por ejemplo, una única central provee más 
del 10% de la demanda. 

La demanda máxima diaria de potencia en Uruguay ronda los 
1000MW. La mayor central térmica tiene una potencia de 125MW, 
por lo que el sistema es muy pequeño y tenemos que tener en 
cuenta el eventual gap. Por otro lado, al ser el sistema tan pequeño, 
la maldición de la dimensionalidad se torna manejable, sobre todo 
teniendo en cuenta la gran capacidad de cálculo con la que se cuenta 
hoy en día. Incluso técnicas de cálculo en paralelo pueden mejorar 
todavía más esta ventaja. 

El presente trabajo analiza condiciones necesarias y suficientes 
para la estabilidad de un algoritmo de optimización dinámica es- 
tocástica para el cálculo del costo del agua en un sistema con una 
única central hidroeléctrica. Se deduce una relación entre el paso de 
integración temporal, el paso de discretización espacial (la altura del 
embalse en este caso) y los caudales máximos de entrada y salida. 
Para mostrar ventajas y desventajas de este enfoque, se presenta una 
simulación del sistema uruguayo con cuatro centrales hidroeléctricas. 
Se verifica que las restricciones en los pasos de discretización son 
muy grandes, lo cual resulta en grandes tiempos de simulación. 

El artículo está organizado de la siguiente manera: En la Sección 
II se presenta el problema de programación dinámica, la deducción 
de la función de costo futuro utilizando una aproximación lineal y 
las distintas formas de implementar las derivadas involucradas. En 
la Sección II, se realiza un análisis de estabilidad asintótica para 
un sistema con un solo embalse y se establece una condición de 
convergencia. Un ejemplo que ilustra el caso en el que hay más de 
un embalse es desarrollado en la Sección IV. Finalmente, el artículo 
se cierra con la exposición de algunas conclusiones y posibles líneas 
de trabajo futuro. 


II. EL PROBLEMA 


Como ya se ha hemos mencionado, se quiere optimizar la 
operación de un sistema eléctrico con una planta de generación 
hidroeléctrica. Introducimos un paso temporal de discretización para 
realizar la simulación. Denotemos por k al instante actual. El estado 
del sistema está dado por el volumen actual V del embalse. Esta 
variable varía desde el nivel mínimo V = O hasta el máximo 
V = Vmaz. Esta variable espacial la discretizamos también, con un 
paso AV; obtenemos así una variable j que varía entre O (embalse 
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vacío) y N (embalse lleno). Se cumple que 


_Vmaz 
N-1 


Construiremos una función de costo futuro FC(V) que nos dice, 
para cada instante de tiempo k y cada nivel posible del embalse j, 
cual es el óptimo costo de operación futura del sistema. 

El Principio de Programación Dinámica [1] nos dice que el costo 
de operación desde el instante k, para un volumen dado V del 
embalse, es igual al costo de operación desde el instante k +1, para 
un volumen V’, más el costo de las decisiones de despacho tomadas 
en el instante k, que llevan el embalse de V a V’. En términos 
matemáticos: 


FC: (V) =B.FC; (V”) + csr (V, us, re) (1) 


AV 





, V =(j-1)AV ,j=1,...,N 


donde 6 representa una tasa de descuento que adapta al presente 
los costos futuros!, V y V’ representan los estados del embalse en 
los instantes k y k + 1 respectivamente y cs; es el costo de las 
decisiones óptimas de operación tomadas en k. Incluye el gasto 
de combustible en las centrales térmicas y el costo asociado a la 
demanda no satisfecha durante la etapa k. Naturalmente, depende 
del estado V del embalse, de las variables de control uz y de ciertas 
entradas aleatorias r, que representan fenómenos tales como los 
diversos caudales de entrada al embalse. Asumimos también que 
conocemos la relación entre el actual nivel V y el nivel final V’. 


£ 
V = f (V, uk, rr) 
Realizamos una aproximación a primer orden de F'C%: 


ƏFCk+ı 
t-a 


Por otro lado, la variación total del volumen del embalse, V’ — V, 
puede escribirse como 


V' -V = Q(V).AT (3) 


FC (V^) = FOr (V) VV -vV 0) 


donde Q; (V) es el caudal neto que ingresa y AT es la duración de 
la etapa. El caudal neto tiene diversas componentes: 


Qk(V)= Qe Qs = Qs — (Qr+Qv + QP) 


QuE denota el caudal entrante e incluye el agua que viene desde 
por el río y el agua de lluvia ;Qr representa el caudal turbinado 
por la central, es decir, agua utilizada para generar energía eléctrica; 
Qv es el caudal vertido y no utilizado para generación; Qp es el 
agua perdida por evaporación y filtración. Combinando (1), (2) y (3), 
obtenemos 


FC4(j) = 6. [Fo0) 





ƏFCk+1 
t-a 


donde hemos simplificado la notación del volumen. A continuación, 
construimos el vector de costos futuros al comienzo de la etapa k, 


cr = [FC (1), FC (2), .. ., FCk(N)] 


(i)-Qk-AT| + csk(j) (4) 


y el vector de costos de operación (usualmente costos asociados a 
las centrales térmicas) en la etapa k: 


Yk F [esk (1), csk(2), ...) csk(N)]” 
Denotemos por q; la matriz diagonal de elementos 


gli i) =Q) , i=1,...,N 


lPor ejemplo, en una base diaria y para una tasa de descuento anual del 


12%, tenemos que 
ab. 


p=( 1 Ja 
—X1+0.12 


Tasas de descuento realistas no deberían exceder el 12%. 


; a əFC l 
Entonces, podemos aproximar la derivada == (j) como com- 
binación lineal de los elementos de Ck+1, dividiendo por el paso 


espacial AV: 





ORC yy 


avy O) = gay Pam 
La ecuación (4) puede escribirse como 
AT 
Cx = 6. (G + Sy hD) .Ck+1 + Yk (5) 


donde D es una matriz que aproxima la derivada. Podemos utilizar 
diversos tipos de matrices D. Por ejemplo, si utilizamos una aprox- 
imación basada en un cociente incremental hacia atrás: 


FCk+1()—-FCk41(i—1) 








LC AV NA 
ƏV FC2(j)-FCpk+41(1) ; 
AA ; j=l 
obtenemos la matriz 
=-1 1 
=l 1 
=1 1 
Didec = i 4 (6) 
=1 1 
—1 1 
E 1 
La repetición de la primera fila se realiza para resolver la condición 
de borde en j = 1. De manera similar, utilizando un cociente 
incremental hacia adelante, resulta la matriz 
—1 1 
-1 1 
Dinc = 2 1 (7) 
=1 1 
=h 
—1 1 


La matriz Daec es adecuada para tratar con el vector de caudales 
salientes qs, dado que cuando el embalse está vacío, no puede salir 
más agua. De manera similar, usaremos Dinc para lidiar con el vector 
qe de caudal entrante. Para una mejor aproximación de la derivada, 
utilizamos una combinación lineal de las matrices Dye. and Dinc. 
La aproximación de primer orden que finalmente utilizamos es: 


AT 

ck = b. |1 + —. 

8 +7 

donde qex y qsx denotan matrices diagonales, constrida a partir de 
los respectivos vectores, similar a la ya mencionada qx. 


(q€x.Dine == qSk-Daec) -Ck+1 + Uk (8) 


III. ANÁLISIS ASINTÓTICO 


En esta Sección, trataremos de encontrar condiciones para la 
convergencia del algoritmo (8). Antes que nada, debemos inicializar 
en k = +00 (o tiempo grande). Por comodidad, revertimos el tiempo 
y consideramos el infinito como k = 0. Consideremos la matriz 


AT 
Aj =P E + AV (dex Dinc — ask Daec)| , Dr = Yk 
Entonces, el algoritmo se expresa como 
Ck+1 = AxCx + bk (9) 


Debemos enfatizar que, en cada etapa, Az y bx denotan las decisiones 
óptimas de despacho, que implican turbinar o vertir caudal del 
embalse. Suponemos conocidos los valores (Q Emaz -máximo caudal 
que puede entrar al embalse- y QSmaz -máximo caudal que puede 
salir del embalse. Denotemos por A y B los conjuntos dentro de los 
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cuales elegimos las matrices Az y by respectivamente. Observemos 
que estos conjuntos son convexos y acotados. 

Consideremos dos condiciones iniciales diferentes co and co, 
es decir, dos valores finales distintos. El algoritmo (9) produce 
dos secuencias {ck} y {c}. Probaremos que, en cada paso, las 
secuencias se aproximan estrictamente. Para medir esta aproximación, 
utilizaremos la norma infinita para vectores: 


lolle = 1 


En el tiempo k + 1, tenemos que 


Ck+1 =  AxCk + bk 
(10) 
Grp = + b, 


La optimalidad de las parejas (Ax, bx) y ( 
consecuencia 


4» Dj.) tiene la siguiente 





Ck+ = As Ck + dr. 








ArCk + bk 


ch = 


(la desigualdad hay que entenderla componente a componente). 
Entonces 


lck+1 = kalloo = |(cr+1); — (c%+1),| 


y ocurre una de las siguientes dos posibilidades 


[(er41); = (ckt1); 
[(Cr+1); a (cka); 


Sabemos que 


N 
[Las (e —d)],|=[ Aey (0:44), 
j=1 


Por lo tanto 





< [Ak (ck — ch)l; 





(11) 











IA 


[As (cx = Ch), 





[Lan (æ -4)],| < 


donde ||A||ı representa la norma matricial inducida por la norma de 
vectores que estamos usando [5]: 


N 
(Al = max, ([Avll=o) = , max aZ 
j=1 


Para cada k, tenemos la desigualdad 





lcr+1 — Ch+1/loo < max (|| 44I|1, 4/11) llcx — Chlloo (12) 


Concluimos que si toda matrizAz € A has [Axl] < 1, las dos 
secuencias generadas por el algoritmo (9) se aproximan a medida 
que el tiempo transcurre (hacia atrás). Recordemos que 


AT 
Ax = b. E + AV (qer Dinc F dsr Dace)| 
Entonces, si consideramos los elementos no nulos de las filas de Az, 
tenemos para la primera fila que: 


AT 


Ar(1,1)=8. [1 y ED - QS:0)]] 


Ar(1,2) = 8.27 [QBw(1) — QS4(2)] 
Para las siguientes filas, 
Arli i1) = p. QSC) 


AT 


Aid =A. fi - 57 [QEx(i) QRO] 


AT 


Arli i +1) = b. zy REO 


La última fila de Az cumple que 


Ax(N,N — 1) = 6.57 [QS:(N) — QEe(N)] 


A(N,N)= 6. [1 Ey [QSA(N) — QEAN) 


Un detalle importante es que la suma de los elementos de las filas de 
Az es siempre igual a 8. Esto implica, en particular, que Az tiene 
un autovector con todas sus componentes iguales. En segundo lugar, 
si elegimos 47 suficientemente pequeño, todos los elementos de Az, 


AV 
son no negativos. Debe ser 
AT QE + QSA(i)] <1 e A 
AT EM) — QSD] <1 (14) 
AT [QSE(N) — QEL(N)] < 1 (15) 


Los únicos términos que pueden presentar problema son Ax(1,2) 
y Ax(N, N — 1). Concentrémonos en el primero. Corresponde a 
la situación del embalse a su mínimo nivel. Entonces, debe ser 
QSx(1) < QEx(1), dado que es una restricción impuesta al 
despacho óptimo (no puede salir más agua que la que entra, pues 
el embalse está vacío). De manera similar, con el embalse lleno, no 
podemos acumular más agua, por lo que el despacho óptimo garantiza 
QSk(N) > QEr(N). Entonces, el cálculo directo da 


Arlı = 8 


y resulta que toda tasa de descuento no nula garantiza la convergencia. 
Podemos también concluir que todos los autovalores de la matriz Az 
están dentro del círculo unitario, ya que [5] 


maxf|Aa, 1) < [Ar ll1 


Hemos probado entonces que la distancia entre dos secuencias 
distintas generadas por el algoritmo (9) tiende a cero en la medida 
que transcurre el tiempo (hacia atrás). Cada una de estas secuencias 
es monótona, ya que representan el costo óptimo de operación, y este 
aumenta al agregar una nueva etapa. Sean a y Kg cotas para los 
elementos de A y B respectivamente (como ya vimos, podemos elegir 
Ka < 1). Consideremos una condición inicial co para el algoritmo. 
En la primera iteración, tenemos que 


cı = Aoco + bo 
De donde 
lIcolloo < |lAoll1-llcollo + llbollo < Kallcollo + Kg 


Para la k-ésima iteración, se cumple que 


k—1 
llexllo < Ilcoll + Xs. (S xs) 


1=0 


Kg 


œ £ [ee] AP 
lczlloo < lleol + 1p 


Entonces, siempre tenemos secuencias monótonas y acotadas, por lo 
que son convergentes. Concluimos pues que para cualquier condición 
inicial, el algoritmo converge a una única solución. 
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TABLE I 


DATOS DEL EJEMPLO DE LA SECCIÓN IV.(**): ALTURA RESPECTO DEL NIVEL 


DEL MAR; (**): LOS VALORES CORRESPONDE AL 50% DE LA PLANTA, DE 


PROPIEDAD URUGUAYA. 






































Bonete Baygorria | Palmar | Salto-UY 
Mínima cota del lago [m] 70 53 36 30 * 
Máxima cota del lago [m] sl 56 44 35.5 * 
Cota de descarga [m] Baygorria Palmar 7.5 5 la 
Capacidad del lago [Hm?] 8210 216 2575 3058 Jia 
Caudal medio de entrada [mi / s] 567 0 290 2358 
Máximo caudal de descarga [m3 /s] 680 828 1373 4200 ER 
Potencia instalada [MW] 155 108 333 945 bl 
TABLE IU 
PARÁMETROS DE OPTIMIZACIÓN DEL EJEMPLO DE LA SECCIÓN IV. 
Bonete | Baygorria | Palmar | Salto-UY 
Máximo caudal de entrada [m3 /s] 567 680 970 2358 
Máximo caudal de salida [m3 /s] 680 828 1373 4200 
Tiempo de llenado del lago (TFL) [días] 168 4 31 15 
Tiempo de vaciado del lago (TEL) [días] 139.7 3.0 21:7 8.4 
Pasos de discretización [n] 10 5 5 5 
Pasos de volumen (V/(n — 1)) [Hm] 912.2 54.0 643.8 764.5 
Tiempo de llenado de un paso de volumen (TFVS) [horas] 446.9 22,1 184.3 90.1 
Tiempo de vaciado de un paso de volumen (TFVS) [horas] | 372.64 18.12 130.26 50.56 























value of water of Baygorria at state x3 (alone) 
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Fig. 1. Valor del agua en Baygorria en el estado x3 - una única central 


hidroeléctrica. (El tiempo está revertido y medido en horas). 


IV. EJEMPLO 


El sistema eléctrico uruguayo tiene cuatro centrales hidroeléctricas 
de generación: ”Bonete”,”Baygorria” y ”Palmar” sobre el río Negro 
y la central bi-nacional de ”Salto Grande” sobre el río Uruguay, 
compartida con Argentina. Las tres plantas sobre el río Negro están 
encadenadas, con Bonete aguas arriba, Baygorria en el medio y 
Palmar aguas abajo. Los parámetros más relevantes de las plantas 
se muestran en la Tabla I. 

A efectos de tener una medida de la estabilidad del algoritmo, 
hemos realizado la optimización con diferentes pasos de tiempo. 
En la optimización, consideramos valores medios para los caudales 
de entrada a los embalses y para las potencias de las distintas 
centrales, es decir que no consideramos efectos estocásticos. Para 
la discretización espacial, tomamos cinco puntos para los lagos de 
Baygorria, Palmar y Salto Grande y diez puntos para el de Bonete. 

Conociendo los valores para los caudales máximos de entrada y 
salida de los lagos, podemos calcular el tiempo de llenado del lago 
(TEL), el tiempo de vaciado del lago (TEL), el tiempo de vaciado 
de un elemento de volumen (TEVS) y el tiempo de llenado de un 


value of water of Baygorria at state x3 (with the others lakes in the system) 
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Fig. 2. Valor del agua en Baygorria en el estado x3 - con todos los lagos 
del sistema uruguayo. (El tiempo está revertido y medido en horas). 


elemento de volumen (TFVL). Estos valores se muestran en la Tabla II 
para los cuatro lagos de Uruguay. El criterio de estabilidad (13)-(15) 
impone un paso temporal menor que el mínimo TFVS y menor que 
el mínimo TEVS. En el ejemplo, el TEVS de Baygorria determina 
que el paso temporal debe ser menor que 18 horas. 

En primer lugar consideremos sólo presente el lago de Baygorria, 
con 3 pasos de discretización. En cada etapa de optimización, el costo 
futuro se calcula para cada uno de los posibles estados del lago x1, 
12, £3, £4 y 5. La Figura 1 muestra los resultados obtenidos para 
el estado x3 para distintos pasos temporales. Para pasos temporales 
menores que 6 horas, los resultados son prácticamente idénticos. 
Para un paso de 12 horas, se aprecian ciertas desviaciones, si bien 
el algoritmo es estable. Para un paso de 18 horas, el algoritmo es 
inestable y los cálculos divergen. 

Cuando consideramos más de un lago en el sistema, es esperable 
que la restricción del paso temporal que calculamos anteriormente ya 
no sea aplicable. La Figura 2 muestra la misma variable que la Figura 
1, pero considerando la presencia de todos los lagos. Observamos aquí 
que las trayectorias correspondientes a pasos temporales de 6 horas 
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o menos son estables y con resultados similares, en tanto que para 
un paso de 9 horas se observa un error importante en el cálculo. 


V. CONCLUSIONES 


En este trabajo hemos deducido un criterio de estabilidad para la 
convergencia de un algoritmo de optimización dinámica estocástica 
para simular el sistema eléctrico y definir así una política óptima 
de despacho. El criterio muestra que hay un compromiso entre el 
paso de discretización temporal, el paso de discretización espacial 
en los embalses y los flujos máximos de entrada y salida de caudal 
de los mismos. En virtud de las simulaciones realizadas, parece ser 
razonable elegir un paso temporal del orden de 5 veces menor que 
la cota teórica, ya que esto parece asegurar un buen control del error 
numérico, además de la estabilidad del algoritmo. Esto es importante 
ya que el algoritmo utiliza una aproximación lineal del problema. 
La cota teórica fue deducida para un sistema con un solo embalse 
y no se aplica directamente a sistemas con más embalses. A tales 
efectos, hemos incluido un ejemplo en el que se considera el sistema 
uruguayo, primero con solo embalse en Baygorria y luego con los 


cuatro embalses que realmente existen. Es notoria la reducción nece- 
saria en el paso temporal para obtener resultados razonables. En la 
actual implementación del simulador, la condición de estabilidad fue 
mejorada usando técnicas de predicción-corrección, para mantener 
acotado el error numérico acumulado. El simulador está disponible 
en forma gratuita en el sitio: http://iie.fing.edu.uy/simsee/. 
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